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^D ■ We propose a novel approach to the Graphene system using a local field theory of 4 dimensional 

(N 

r^ . QED model coupled to 2+1 dimensional Dirac fermions, whose velocity is much smaller than the 

speed of light. Performing hybrid Monte Carlo simulations of this model on the lattice, we compute 

00 ■ 

the chiral condensate and its susceptibility with different coupling constant, velocity parameter and 

"j^ . flavor number. We find that the chiral symmetry is dynamically broken in the small velocity regime 

o ' and obtain a qualitatively consistent behavior with the prediction from Schwinger-Dyson equations. 
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I. INTRODUCTION 

The Efraphene material which is composed of two-dimensional carbon atomic seat has 

n 

attracted much attention from its specific property [1|]. The high conductivity of Graphene 
has understood by the analogical picture to "relativistic" massless Dirac fermion (quasipar- 
ticle) bounded in 2+1 dimension which has been derived in based on tight-binding (TB) 
model at low energy 2, |3|. In the effective model with massless Dirac fermion, compared 
with quantum electromagnetic dynamics (QED) in particlephysics, the velocity of fermion 
is 2-orders magnitude smaller than the speed of light (2|, |3[. Recent experiments suggest 



fl 



that this magnitude may be changed by a factor [J] due to significant quantum corrections. 
The spontaneous gap generation at semimetal-insulator transition point (Mott transition) 



of Graphene has been predicted by Schwinger- Dyson equation and large N model [SMlO] 
in the presence of Coulomb interaction between Dirac fermion before Graphene confirmed 
in experiments so far, however experimentally there has been no conclusive evidence that 
semimetal-insulator transition of suspended monolayer Graphene occurs in low temperature 



even in 20 K [ll| (as opposed to this bilayer Graphene posses energy gap structure [12, Il3j). 
It may be that the model is missing an important feature of the dynamics or that the approx- 
imations used for the theoretical prediction lead a wrong conclusion. In order to clarify this 
point it is necessary to carry out an exact calculation based on a solid theoretical approach. 
The above effective Dirac fermion system possesses chiral symmetry which could be spon- 
taneously broken by strong effective Coulomb interaction due to enlarged effective fine struc- 
ture constant Og = e^ /{At:vf) ~ 2.2 at small velocity vp ~ 1/300 estimated as vp = 3ta/2 [2(] 
in the natural unit. Here, t denotes leading hopping parameter of TB model and a denotes 
distance among the nearest-neighboring carbon atoms. The semimetal-insulator transition 
of Graphene is considered as a consequence of the energy gap generation in the phase of 
chiral symmetry breaking (xSB). This is analogous with the mass-gap generation in 2+1 
dimension QED system suggested in 1/N expansion [1^, Schwinger-Dyson equation 15| 20 
years ago and also several attempts with Monte-Carlo simulation in non-Compact lattice 



QED |16l-ll8|. However, the present system possesses crucially different points from the 2+1 
dimension QED system. Firstly the quasiparticle has different velocity from speed of light in 
spatial 2 dimension brane and secondly there is Coulomb interaction in spatial 3 dimension 
but not in spatial 2 dimension. The condensate for scalar bilinear operator {ipip) is also as 



an order parameter of xSB at low temperature in this QED system. 

Note that "chiral" symmetry in 2+1 dimension system is defined for the four component 



Dirac spinor. According to reference [19|, four component Dirac spinor is conventionally 
described as ^ = {■ip^~^,ip^~^,ip^~ ,ip^~) where ip^'^"^ denotes electron wave function (Bloch 
state) located on sublattices (A,B) in hexagonal Graphene lattice and two independent 
energy grand points (valleys) (±). The subscript means the spin indices of electron. In 
this notation the chirality of quasiparticle corresponds to valley indices which can be dis- 
tinguished by 75 projection as well as the particle physics. The gamma matrix is given by 
tensor structure as 70 = o-o^hx2, li = —i(T2®0'i {i = 1, 2, 3) and 75 = «7o7i7273 constructed 
by spin and sublattice-valley indices. In the case of monolayer Graphene, chiral symmetry 
can be regarded as global "fiavor" UjA ) sy mmetry of quasiparticle whose generator is given 



by {1,75, Z73, [73,75]/2} ® crj=o,i,2,3/2 [19[ in which the first group generates the rotation 
operator for sublattice and valley points and second one generates spin rotation (in general 
multilayer Graphene is described in U(2Nf) symmetry by taking account of layer number 
A'laycr as fiavor Nf = N^pin x A^'iaycr assuming there is no interaction between layers). Using 
the electron wave function the chiral condensate can be expanded into 



{^^p) = Y, I [v^^+ V" + ^^" V"" + /i-c] (1) 



which is indeed mixture of different chirality of Graphene [19| . 

In the leading order of extended TB model into "non-relativistic" QED including Coulomb 
interaction between quasiparticle, there have been many model calculations of chiral con- 
densate near critical point. For example, the Schwinger- Dyson equation predicts that there 
is a critical coupling Ur above which the system has a nonzero condensate. In the mean-field 

n 

approximation, ac is predicted to be Oc = 0.5 [20||. When one-loop vacuum polarization ef- 
fect is taken into account, the gap equation provides stronger critical coupling Uc = 0.92 |8|. 
The similar critical value has been also obtained by lattice calculation in the Monte-Carlo 



simulation 
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or strong coupling expansion [23| using "non-relativistic" QED action. 
In this paper we propose a new strategy for non-perturbative study of the semimetal- 
insulator transition of Graphene with QED model. The important difference from the frame- 
work of "non-relativistic" QED is that we consider QED action manifestly including both 
gauge invariant interaction and velocity contribution in the framework of local field theory. 
In our realistic setup it is possible to rigorously investigate the velocity contribution to xSB 



phenomena and compare to actual Graphene. We numerically show that critical behavior 
when changing velocity and flavor number is qualitatively consistent with Schwinger-Dyson 
equation. 

This paper is organized as the follows: In section [IT] we explain the model of relativistic 
QED action with fermion brane and in section IIIII we show the detail of setup of our lattice 
action and simulation. In section IIVI we show the behavior of chiral condensate and sus- 
ceptibility when changing the parameters of coupling constant, fermion mass, velocity and 
flavor number. In section |V] we summarize and discuss the future works. 

II. MODELING THE RELATIVISTIC QED 

In the graphene system, since the velocity of the effective Dirac fermion is extremely small, 
the following action with only the instantaneous Coulomb interaction has been adopted |24l |: 



•S'couiomb = dtd xijj 



idtjt + ivF{dxlx + idy'^y) 



^ 



e^^ft - 1'\ - 

dtd'^xdt'(fx'^p-ftijit,x) — ^ j^jt^(t,x'). (2) 

87r|a; — x'\ 

Introducing the scalar potential as an auxiliary fleld, one could equivalently describe the 

action as 

1 



Snr = 7^-^ / dtd^xE"^ + / dtd'^x ^ 



i{dt + i(f))'yt + ivF{dxlx + dy'yy) ^p, (3) 



where E = — V0. 

QED action with fermion bounded on 2+1 dimensional "brane" under gauge invariance 
in the continuum theory is straightforwardly written as 

S = - I dtd^x{E^ + B^)+ I dtd^xilj[iDt-ft + iv{D^-fy + Dy-fy)]-^ (4) 

with fermi- velocity t" as a coefficient of gauge covariant derivative Di = di + iAi and coupling 
constant /3 = 1/e^ (here we distinguish velocity parameter "w" from renormalized (physical) 
velocity "wi?" as explained below.). Electric and magnetic fleld can be described as Ei = Fit, 
Bi = eijkFjk/2 {i = X, y, z) with fleld strength F^^, = d^A^ — d^A^, which is deflned on 3+1 
dimension, while fermion fleld ip is bounded on 2+1 dimension. Second using simultaneous 
scale transformation for temporal coordinate and temporal gauge fleld normalization, 

t -)■ t/f , Aq -)■ Aqv, (5) 



Eq.dl]) is represented as 



s = ^ 

2 



j,uM^E^,.-^B').JMS.ij:iO,,,,, 



dt<fx{vE'^ + v"^B^)+ I dtSx%l) Y^ iD^-ff,ij, (6) 

thus fermi-velocity shows up in the anisotropy of gauge couphng. As a consequence of 
rescahng we can equivalently deal with kinematic term of massless fermion field as the 
standard QED. Adopting the rescaled photon field for spatial direction Ai -^ AiV, Eq.(l6]) 
can be described as 



S 



vf3 

T 



dtd^x \Y,i^^Ao - vdoA.y Yl +(^^^i - ^^-^^ 



^,J 



dtd^xij\i{do + iAo)-ft + ^ i{di + wAi)-fi 



^, 



(7) 



i=x,y 



Taking the non-relativistic limit, which is interpreted as limit of w — )■ while v(3 is fixed, 
the spatial photon field can be decoupled in the path integral. Dropping the spatial photon 
term, non-relativistic action can be derived as 



^1 



NR 






dtd^xE'^ + / dtd'^x ^ 



iDtjt + id^-f^ + idy-fy 



^, 



(8) 



so that in this limit time- rescaled version (t — )■ t/v and — )■ v(j)) of the action in Eq.(l3]) is 
reproduced. 

Here we notice that the relativistic action ([6]) is a renormalizable local quantum theory. 
Therefore, the scaling of velocity can be rigorously treated without the contamination from 
irrelevant operator violating gauge symmetry. The bare fermi velocity v has a renormaliza- 
tion effect due to quantum correction to ultraviolet divergence. Indeed the renormalization 



25l 



26| (and also there is similar dis- 



group study from relativistic perturbative calculation 

cussion in large N 27| and renormalization group 28l-l30| in spite of non-relativistic action) 
pointed out that the velocity parameter behaves a logarithmic-divergent scaling due to vir- 
tual fermionic loop correction. It turns out that the the renormalized fermi velocity vp has 
an infrared fixed point at vp = 1 due to the restoration of Lorentz symmetry. Perturbative 
analysis has also pointed out that there is an infrared unstable (ultraviolet) fixed point of 
the running of fermi-velocity above the experimental one (see Figured]). However, as no- 



ticed in 
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, l26| , such a fixed point is in the strong coupling region where the perturbation 



theory is not reliable. In addition that recent experimental study using effective spectrum of 
electron-hole pair in suspended Graphene on Si02 wafer [^ also shows the effective velocity 




FIG. 1: The sketch of renormahzation flow of renorniahzed fermi- velocity vp and effective coupling 
a with avp =const. In the perturbation theory there is unstable infrared fixed point at Ofix which 
is shown as straight dashed-line. Experimental point represented as cross symbol at vp ~ O(10~^). 
Our lattice calculation aims to prove the wide region of phase structure {a,VF) shown as the gray 
zone beyond the region perturbatively reliable. 

of Graphene around a few particle region has qualitatively similar behavior of renormalized 
velocity in the perturbation theory. Therefore, investigation of running behavior of fermi- 
velocity using non-perturbative calculation with relativistic action would be very important 
task towards more rigorous discussions than renormahzation group study. 

III. LATTICE SIMULATION OF THE RELATIVISTIC QED 

In our Monte-Carlo study, we employ the conventional non-compact gauge action 

4D 



/3, 



S9 = Y1 W^W.ix) - 9A(x))' + -{d^e,{x) - dM^)) 



(9) 



with anisotropic gauge coupling and vector field 6^. This summation is performed in 4 
dimensional coordinate space-time. In this action there are two input parameters, /3 and v, 
which correspond to the bare gauge coupling constant and velocity. For the fermion action 
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on the lattice, we use the staggered type 

3D 

5/ = ^ |^^?7^(x)|x(a;)t/^(x)x(a; + /i) - U\x - fi)x{x - fi)j + mx{x)x{x) (10) 

X fj, 

where hnk variable is defined as U^{x) = exp(z0^(x)), and rj^^x) = n^^~^(— 1)* is Kawamoto- 
Smit phase factor. In the above equation the staggered fermion is defined in 3 dimensional 
field, and thus the summation is defined in 3 dimensional coordinate space-time. For com- 
parison the size of z-axis, which is perpendicular coordinate for fermion brane, is set to 
same L^ = 8 as 21|. Boundary condition of fermion in temporal direction is anti-periodic 
boundary associated with finite temperature, and the other directions are set to periodic 
boundary. In this setup gauge field propagates along z-axis without interaction of fermion. 
We insert the small mass term with parameter ttt, as a probe to investigate chiral dynamics 
near chiral limit. 

Conveniently staggered fermion in 3 dimensional coordinate automatically satisfies a same 
fiavor U{A) symmetry as monolayer Graphene without additional counter term in the con- 



tinuum limit 2l| (due to the lattice artifact it is broken to U{1) ^ U{1)). Lattice size is 



30^ X 20, and we perform the Hybrid Monte Carlo (HMC) method with Omelyan integrator 



using parameter Ac = 0.1931833 [311 and Hasenbusch accelerator |32[ with mass parameter 
m = 0.05 to generate gauge configurations with dynamical fermion. We tune the Leapfrog 
time-stepping parameter 6t = l/Nt, where Nt is number of step in unit HMC trajectory, to 
be that acceptance rate is satisfied with more than 60% in our simulation. From the practical 
point of view, the rescaling (|5]) has also an advantage of simultaneously covering the region 
of strong electric coupling and low temperature. Since temporal size is 1/v times larger than 
original action (jl]), the lattice simulation curries out naively equivalent to be at 1/v times 
smaller temperature 33|] with same computational cost. Note that the physical quantities 
obtained in this action depends on the effective strong coupling a = 1/{4ttPv) rather than 
gauge coupling constant (3 since the quantum correction from vertex with a is dominant 
contribution in this system if velocity is enough small. Correspondingly the loop correction 
to gauge coupling constant is higher order effect while velocity has logarithmic divergence 
(f (/i) = f (1 + a/41n(A/u) with cut-off A) from at least one-loop correction according to the 



perturbative analysis [25[. It turns out that observable are mostly controlled by the theo- 
retical parameter v instead of gauge coupling. Furthermore we introduce the Dirac "mass" 
term mi/jilj which preserves discrete symmetry (Parity, time- reversal and charge conjugate) 



but not chiral symmetry. Scalar condensate of quasiparticle {^ip) corresponds to order pa- 
rameter of chiral transition, and associated from Goldstone theorem the massless Goldstone 
mode (pion) appearing into QED model plays a significant role of the gap generation in 



Graphene system 34 1 



IV. RESULT OF LATTICE SIMULATION 

Here we investigate the chiral breaking phenomena with chiral condensate a = (xx) ^^^ 
its susceptibility Xm = da /dm behaved as a function of /3f , v and m. The disconnected 
diagram appearing in chiral condensate a and chiral susceptibility Xm is calculated by the 
noise method with 200 noise vector. We sample configurations at every 20 HMC trajectories, 
and employ the Jackknife error estimate with 10 bin-size. Total statistics used in this analysis 
is around 600-1500 configurations, especially more than 1000 configurations in the strong 
coupling region (/3t> < 0.06). 

Chiral condensate is an order parameter of spontaneous chiral symmetry breaking in 
which the discontinuous point along the coupling constant regards as second-order critical 
point in the limit of massless on infinite volume. At the finite mass and lattice volume, if 
there is critical point associated with spontaneous chiral symmetry breaking in this system, 
chiral condensate above critical coupling a^ remains finite value, otherwise below ac chiral 
condensate advances toward zero when approaching ?7i = 0. Figure [2] clearly shows the 
expected critical behavior that, when mass parameter is decrease, a at weak coupling region, 
which is at the large /3w, goes to zero while at strong coupling region a remains finite value, 
and critical coupling is found to be Pv = 0.05-0.06 at fixed v = 0.1. Chiral susceptibility also 
shows significant behavior around such critical point. When m goes down close to zero, peak 
position of Xm grows up and shifts from weak coupling to strong coupling region. We expect 
the singularity for critical phenomena at massless point will appear at I3v =0.05-0.055, which 
corresponds to critical coupling Oc =1.45-1.59. This value is similar order of the prediction 
in g ap equation quoted to ac = 0.92 |8| or in non-relativistic lattice calculation as ac = 1.02 
2l| . Figure [3] also shows a significant change as that below ac we can see that a approaches 
to zero at chiral limit, especially near a^ its non-linear behavior, while above ac u does not 
vanish. At very weak region in which the coupling anisotropy is negligible since /3 is too 
large, a can be described as a perturbative line, and actually in Figure [3] mass dependence 
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FIG. 2: Chiral condensate (left) and chiral susceptibility (right) of quasiparticle as a function Pv 
which corresponds to effective coupling constant a at fixed velocity in w = 0.1. 
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FIG. 3: Mass dependence of chiral condensate at several gauge coupling /3 at fixed velocity in 
V = 0.1. Linear-dashed line shows the leading perturbation result. 



of a approaches to the leading perturbative line obtained in lattice perturbation. 

To see the contribution of fermi velocity and vacuum polarization of bounded fermion 
into critical phenomena, we compare the peak position of chiral susceptibility at t; = 0.1 and 
V = 0.05 (in addition to v = 0.03 for full calculation) with and without dynamical fermion 
at mass parameter m = 0.0025 in Figure IH We notice two points: 

1. First, our lattice calculation shows that critical point in full QED is moved to strong 
region in quenched QED Pv ~ 0.093 corresponding to a^ ~ 0.8, which is about 40-50% 
strongly shift of coupling constant. Our result is qualitatively consistent behavior with 
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FIG. 4: Comparison with and without dynamical quasiparticle at different velocity parameter for 
chiral susceptibility. 

result in gap equation including frequent vacuum polarization effect of 2+1 dimensional 
fermion, quoted to be 46% shift from a^ = 0.5 (quench) |20| to ac = 0.92 



8| 



2. Secondly, decreasing the velocity from v = 0.1 to v = 0.05 we see that the critical 
point a stays almost constant but shifts towards stronger coupling for only about 
5-7% in both cases of full QED; peak position of Xm is moved as f3v ~ 0.062(t; = 
0.1) —7- 0.058(t' = 0.05), and quenched QED; peak position of Xm is moved as f3v ~ 
0.093(f = 0.1) — 7- 0.088(f = 0.05). Linearly extrapolating into v ~ 1/300 using 
ac{v) = ac{v = 0) + cv assuming c is constant, we find that the critical point is 
ac{v = 1/300) ~ 1.5 in the case of full QED {adv = 1/300) ~ 0.9 in quenched QED), 
which is comparable to TB model {a ~ 2.2) |2| rather than quenched one. Assuming 
that the bare velocity v and the renormalized velocity vp is almost the same, our 
results suggest the physical point of the Graphene system exists within the symmetry 
broken phase. To establish the position of the physical point in the phase structure, 
we have to carry out the renormalization of the effective coupling constant as well as 
the fermi velocity, which will be left a future work. 
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V. DISCUSSION 

In this paper we have performed the lattice simulation using 4 dimensional non-compact 
QED with 2+1 dimensional fermion brane including velocity parameter into anisotropic 
gauge coupling between electric and magnetic fields. Our lattice simulation clearly shows 
that a consistent behavior associated with chiral symmetry breaking at critical coupling 
Oc =1.45-1.59 at t" = 0.1, and furthermore from comparison with different velocity param- 
eter we suggest possibility of quantum correction to the critical phenomena depending on 
running of velocity parameter in the infrared region. We assert that in our relativistic QED 
model related to "realistic" Graphene system which includes the fermi-velocity less than 
1% speed of light and electromagnetic interaction among electron, the spontaneous chiral 
symmetry breaking occurs due to small velocity parameter, which plays important role to 
lead to the strong dynamics of this system rather than gauge coupling as suggested by per- 



turbative argument 25, 261. Qualitatively the result obtained in our model is similar to 
model prediction of gap equation [3, |8| . We argue that our model is non-trivial extension of 
QED model under gauge invariance to relativistic form in contrast to non-relativistic case 



211 ]. Although we have not comprehended the explicit relation between our model and 



Graphene in low energy region yet, however surprisingly our model also shows the critical 



m 



m. 



behavior of quasiparticle condensate as well as in the non-relativistic calculation |21|, 

This is strong evidence of that our model realizes the natural extension to relativistic theory 

of Graphene model beyond TB approximation. We are now proceeding to take into account 



renormalization of velocity and continuum limit [38| using more detailed parameter search 
under way. This is a feasible study for one of the most challenging topics of theoretical 
prediction of semimetal-insulator transition point of suspended Graphene. 
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